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KRYLOV APPROXIMATION OF ODES WITH 
POLYNOMIAL PARAMETERIZATION 

ANTTI KOSKELA*, ELIAS JARLEBRING*, AND MICHIEL E. HOCHSTENBACHt 

Abstract. We propose a new numerical method to solve linear ordinary differential 
equations of the type ^ {t, e) = A{e) u{t, e), where A ; C —>■ is a matrix polynomial with 

large and sparse matrix coefficients. The algorithm computes an explicit parameterization 
of approximations of u{t, e) such that approximations for many different values of e and t 
can be obtained with a very small additional computational effort. The derivation of the 
algorithm is based on a reformulation of the parameterization as a linear parameter-free 
ordinary differential equation and on approximating the product of the matrix exponential 
and a vector with a Krylov method. The Krylov approximation is generated with Arnoldi’s 
method and the structure of the coefficient matrix turns out to have an independence on 
the truncation parameter so that it can also be interpreted as Arnoldi’s method applied 
to an infinite dimensional matrix. We prove the superlinear convergence of the algorithm 
and provide a posteriori error estimates to be used as termination criteria. The behavior of 
the algorithm is illustrated with examples stemming from spatial discretizations of partial 
differential equations. 

Key words. Krylov methods, Arnoldi’s method, matrix functions, matrix exponential, 
exponential integrators, parameterized ordinary differential equations, Erechet derivatives, 
model order reduction. 
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1. Introduction. Let Aq, Ai, ..., An € be given matrices and con¬ 
sider the parameterized linear time-independent ordinary differential equation 

Ovj 

(1-1) e) = A{e) u{t, e), -(/(O, e) = uq, 

where A is the matrix polynomial A(e) := Aq -|- sAi -!-•••+ An- Although 
most of our results are general, the usefulness of the approach is more explicit 
in a setting where N is not very large and the matrices Aq, ..., An are large 
and sparse, e.g., stemming from a spatial finite-element semi-discretization of 
a parameterized partial-differential equation of evolutionary type. 

We present a new iterative algorithm for the parameterized ODE dLU), 
which gives an explicit parameterization of the solution. This parameterization 
is explicit in the sense that after executing the algorithm we can find a solution 
to the ODE (II. ip for many different values of e and t > 0 without essential 
additional computational effort. Such explicit parameterizations of solutions 
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are useful in various settings, e.g., in parametric model order reduction and in 
the field of uncertainty quantification (with a single model parameter); see the 
discussion of model reduction below and the references in [3]. 

The parameterization of the solution is represented as follows. Let the co¬ 
efficients of the Taylor expansion of the solution with respect to the parameter 
e be denoted by co(f), ci(t), ..., i.e., 

OO 

( 1 . 2 ) u{t,e) = exp{tA{s))uo = E 

e=o 

As exp(tA(e)) is an entire function of a matrix polynomial, the expansion 
dEl exists for all e € C. 

Consider the approximation stemming from the truncation of the Taylor 
series (US!) and a corresponding approximation of the Taylor coefficients 

k-l 

Uk{t,£) := ^e^Q(t) 

e=o 

k-l 

(1.3) PS =■ Ukit, e). 

1=0 

Our approach gives an explicit parameterization with respect t of the approxi¬ 
mate coefficients Co(t),... ,Cfc_i(t) which, via (II.3p . gives an approximate solu¬ 
tion with an explicit parameterization with respect to e and t. 

The derivation of our approach is based on an explicit characterization of 
the time-dependent coefficients Co{t),.. . We prove in Section [2] that 

they are solutions to the linear ordinary differential equation of size nm, 


d 

co{t) 


co{t) 


co(0) 


no 

0 

dt 

Cm—1 {t) 

— 

Cm—1 (0 


Cm—1 (0) 


0 


The matrix Lm in (11.41) is a finite-band block Toeplitz matrix and also a lower 
block triangular matrix. 

Since (jl.4p is a standard linear ODE, we can in principle apply any numer¬ 
ical method to compute the solution which results in approximate coefficients 
Co(t),. .. ,Cm-i(f). Exponential integrators combined with Krylov approxima¬ 
tion of matrix functions have recently turned out to be an efficient class of 
methods for large-scale (semi)linear ODEs arising from PDEs [11[I12| . See also 
m for a recent summary of exponential integrators. Krylov approximations of 
matrix functions have a feature which is suitable in our setting: after one run 
they give parameterized approximations with respect to the time-parameter. 

Our derivation is based on approximating the solution of (II.4p . i.e., a prod¬ 
uct of the matrix exponential and a vector, using a Krylov method. This is 
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done by exploiting the structure of the coefficient matrix Lm- We show that 
when we apply Arnoldi’s method to construct a Krylov subspace corresponding 
to (HID, the block Toeplitz and lower block triangular property of result 
in a particular structure in the basis matrix given by Arnoldi’s method. 

The structure of is such that, in a certain sense, the algorithm can be 
equivalently extended to infinity. For example when N = 1, the basis matrix 
is extended with one block row as well as a a block column in every iteration. 
This is analogous to the infinite Arnoldi method which has been developed for 
nonlinear eigenvalue problems |15j and linear inhomogeneous ODEs jl6) . This 
feature implies that the algorithm does not require an a priori choice of the 
truncation parameter m. 

We prove convergence of the algorithm (in Section [3D and also provide a 
termination criteria by giving a posteriori error estimates in Section 01 

The results can be interpreted and related to other approaches from a 
number of different perspectives. From one viewpoint, our result is related 
to recent work on computations and theory for Frechet derivatives of matrix 
functions, e.g., uni da 08]. As an illustration of a relation, consider the special 
case N = 1. The first-order expansion of the matrix exponential in (|1.2I) and 
m Chapter 3.1] gives 

u{t, e) = exp(f(Ao + eAi)) uq 

= exp{tAo)uo + Le^p{tAo,etAi)uo + o(|e| |t| ||Ai||), 

where Lexp is the Frechet derivative of the matrix exponential. Since the 
Frechet derivative is linear in the second parameter, the first coefficient is ex¬ 
plicitly given by ci{t) = Lexp(tAo, tAi) uq- The higher order terms C 2 ,C 3 ,... 
have corresponding relationships with the higher order Frechet derivatives. An 
analysis of higher order Frechet derivatives is given in [20] . In contrast to 
the current Frechet derivative approaches, our approach is an iterative Krylov 
method with a focus on large and sparse matrices and a specihc starting vec¬ 
tor, which unfortunately does not appear to be easily constructed within the 
Frechet derivative framework. 

The general approach to compute parameterized solutions to parameter¬ 
ized problems is very common in the field of model order reduction (MOR). See 
the recent survey papers EllSI. In the terminology of MOR, our approach can 
be interpreted as a time-domain model order reduction technique for parame¬ 
terized linear dynamical systems, without input or output. Parametric MOR 
is summarized in [3|; see also ini[23]. Our approach is a Krylov method to 
compute a moment matching approximation in the model parameter £. There 
are time-domain Krylov methods, e.g., those described in PhD thesis [6]. To 
our knowledge, none of these methods can be interpreted as exponential inte¬ 
grators. 

We use the following notation in this paper. We let denote the num¬ 
ber of elements in the set S, and vec{B) denote vectorization, i.e., vec(i?) = 

3 


[bj,..., G C"'*', where B = [bi,...,bk] G By we indicate the 

identity matrix of dimension n. The set of eigenvalues of a matrix A is de¬ 
noted by A(74) and the positive integers by N+. The logarithmic norm (or 
numerical abscissa) —)■ M is defined by 

(1.5) /r(A) := max j^A G M : A G A 

2. Derivation of the algorithm. 

2.1. Representation of the coefficients using the matrix exponen¬ 
tial. To derive the algorithm, we first show that the time-dependent coeffi¬ 
cients Co(t), Cm-i{t) are solutions to a linear time-independent ODE of 
the form (|1.4D . i.e., they are explicitly given by the matrix exponential. 

Theorem 1 (Explicit formula with matrix exponential). The Taylor co¬ 
efficients Co(t),... ,Cm-i(t) in (11.21) are explicitly given by 


A +A* 


(2.1) vec(co(t),... , Cra-i{t)) = exp(tLm) Uq, 
where 


( 2 . 2 ) 


Lm . — 


Aq 

Ai 


A 


N 


A 


N 


Ai Aq 


G 


and uq = 


Uq 

0 


G C" 


and N = min(m — 1, N). 

Proof. The proof is based on explicitly forming an associated ODE. The 
result can also be proven using a similar result |18[ Theorem 4.1]. We give here 
an alternative shorter proof for the case of the matrix exponential, since some 
steps of the proof are needed in other parts of this paper. Differentiating (11.21) 
yields that for any j > 0, 


(2.3) 


1 

jl de^ 


u{t, e) 


e=0 


Cj{t). 


By evaluating (12.31) at t = 0, and noting that m(0,£) = uq is independent of e, 
it follows that co(0) = uq and 0^(0) = 0 for i > 0. The initial value problem 
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(HU) and the expansion of its solution (HU imply that 


1 d 


_ 1 
e=o ~ J! 


N 


u{t,£) 


\i=0 


£=0 


N 


(2.4) 




i=0 


£=0 


£ = 0 


N oo 

EE 

i=0 1=0 


1 

j\ 


S+^ 


£•=0 


AiCi{t) 


min(A^j') 

= ^ AiCj-i{t). 

i=0 

From (I2U and (EU it follows that the vector vec(co(f),..., Cm-i{t)) satisfies 
the linear ODE (11.41) with a solution given by (12.11). □ 


Algorithm 1: Infinite Arnoldi algorithm for polynomial uncertain ODEs 
Input : uq € C"', Aq,. .. ,A]\f G 

Output: Matrices Qp G i))xp and Hp G representing 

approximations of the coefficients cq, ..., Cp_i via (|2.7I) 


1 Let /3 = lluoll, Qi = uo/P, ^0 = [ ] 

for f = 1, 2,... ,p do 


2 

3 

4 


Let X = Q{:,£) G c^+V-^)nN 
Compute y := vec(yi,..., yi+(£_i) 7 v) 


Let 0^ := 





G with (12.6p 


5 

6 

7 

8 

9 

10 


Compute h = Q*w 
Compute yA_'-= y — Qjh 
Repeat Steps [SHS] if necessary 
Compute a = ||y±|| 


Let = 
Let Qf+i : 


Ke- 1 h 
0 a 

= [Q^,y±/a] G c^^+ii-i)nN)xie+i) 


end 

11 Let Hp G be the leading submatrix of H^p G R(p+i)^p 


2.2. Algorithm. Theorem [T] can be used to compute the coefficients C£(t) 
if we can compute the matrix exponential of Lm times the vector uq. We use 
a Krylov approximation which exploits the structure of the problem. See, e.g., 
mm for literature on Krylov approximations of matrix functions. 

The Krylov approximation of v{t) = exp{tB)vo, consists of p steps of the 
Arnoldi iteration for the matrix B initiated with the vector vq. This results in 
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the Arnoldi relation 


BQp — Qp+iH.p') 

where G cip+ilxp is a Hessenberg matrix and Qp is an orthogonal matrix 

spanning the Krylov subspace ]Cp{B,vo) = s]pan{vo, Bvq, ..., The 

Krylov approximation of exp{tB)vQ is given by 

v{t) = eyip{tB)vo « Qpexp(tHp)ei ||uo||, 

where Hp G is the leading submatrix of ^p, and ei is the first unit basis 

vector. 

The only way B appears in the Arnoldi algorithm is in the form of matrix 
vector products. Moreover, the Arnoldi algorithm is initiated with the vector 
vq. Suppose we apply this Arnoldi approximation to (12.111 . In the first step we 
need to compute the matrix vector product 

(2.5) Lm vec{uo, 0,..., 0) = vec(Ao mq, ■ ■ ■, ^o, 0,..., 0), 

which is more generally given as follows. 

Lemma 2 (Matrix vector product). Suppose x = vec(xi ,... ,Xj,0,... ,0) = 
vec{X) G C”™, where xi,... ,Xj G C”' and m > j + N. Then, 

LmX = vec(yi,..., yj+N, 0,... , 0), 


where 

min(Af,£—1) 

(2.6) yi= ^ AiXi-i, ^ = + N. 

z=max(0,£—/c) 


Proof. Suppose S G is the shift matrix S := which 

satisfies (S'*)'^ = We have 

N N 

LmX = ^(5* ® Ai) vec(X) = ^ veciAX{SY) 

2 = 0 2=0 

N m—i 

= EE vec{AiXe£eJ_^_^). 

i=0 £=1 


Note that Xe£ = X£ if £ < j and Xe£ = 0 if £ > j. Hence, by using the assump¬ 
tion m > j+X, and reordering the terms in the sum we find the explicit formula 

A.rr. \ _ s:^j+N ,£-1) a T\ 


LmX = E *=0 ELi vec(AiX£e^+J = ° 

Since the Arnoldi method consists of applying matrix vector products and 
orthogonalizing the new vector against previous vectors, we see from (12.5p that 
the second vector in the Krylov subspace will consist of A^ -f 1 nonzero blocks. 
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Repeated application of Lemma [2] results in a structure where the jth column 
in the basis matrix consists of (j — l)A/^+1 nonzero blocks, under the condition 
that m is sufficiently large, ft is natural to store only the nonzero blocks of 
the basis matrix. By only storing the nonzero blocks, the Arnold! method for 
(j2.ip reduces to Algorithm [TJ 

Note that our construction is equivalent to the Arnold! method and the 
output of the algorithm is a basis matrix and a Hessenberg matrix which to¬ 
gether form the approximation of the coefficients cq , ..., Ck-i 

(2.7) vec(co(t),... ,Cfc_i(t)) « vec(co(t),... ,Cfc_i(t)) := Qpexp(tL/p)ei \\uo\\, 

where by construction k = N(p — 1). The approximation of the solution is 
denoted as (m, i.e.. 


fc-i 

Uk,p{t) ■■= 

e=o 

where we added an index p to stress the dependence on iteration. A feature 
of this construction is that the algorithm does not explicitly dependend on 
m, such that it in a sense can be extended to infinity, i.e., it is equivalent 
to Arnoldi’s method on an infinite dimensional operator. The result can be 
summarized as follows. 

Theorem 3. The following procedures generate identical results. 

(i) p iterations of Algorithm{J\ started with uq and Aq, .. ., A]\f; 

(ii) p iterations of Arnoldi’s method applied to with starting vector ei® 
uq € C"'™' for any m > Np; 

(Hi) p iterations of Arnoldi’s method applied to the infinite matrix L^o with 
the infinite starting vector ei ®Uo & C°°. 

3. A priori convergence theory. To show the validity of our approach 
we now bound the total error after p iterations, which is separated into two 
terms as 

(3 1 ) “ ^N{p-i),p{t, e) II 

< errK,jv(p_i),p(t,e) +errT,Ar(p_i)(t,e), 

where 

(3.2) errK,fc,p(t,e) := \\uk,p{t,e) -Uk{t,e)\\ 

(3.3) errT,fc(t,e) := \\u{t,e) -Uk{t,s)\\. 

A bound of errK,fc,p) which corresponds to the Krylov approximation of the 
expansion coefficients Co,...,Cfc_i, is given in Section l3.ll and a bound on 
errx,A:) which corresponds to the truncation of the series, is given in Section 
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After combining the main results of Section 13.11 and Section 13.21 in particular 
formulas (|3.8I) and (|3.9I) . we reach the conclusion that 


(3.4) errp(t,e) < 

Ci{t,e) 2^ -, n —-lko||+2^ 

i=o 


1 — |e 


i2N{p-l) 


(ta)^e 




{p + £-2)l 


1 — lei 


p\ 


■IKII, 


where a and 7 are given in (13.6p . and C'i(t, e) and C' 2 (t,e) are given in (|3.10p . 
Due to the factorial in the denominator of (|3.4p . for fixed e and t > 0, the total 
error approaches zero superlinearly with respect to the iteration count p. 

3.1. A bound on the Krylov error. We first study the error generated 
by the Arnoldi method to approximate the coefficients co{t),... We 

define 


(3.5) Ek,p{t) = [co{t),...,Ck-i{t)] - [co(t),... ,Cfc-i(t)]. 


where cq, ..., are the approximations given by the Arnoldi method, i.e., by 
the vector 


Ck{t) : = 


Co{t) 




= Qpex.p{tHp)ei ||uo||- 


Using existing bounds for the Arnoldi approximation of the matrix exponen¬ 
tial |B], we get a bound for the error of this approximation, as follows. 

Lemma 4 (Krylov coefficient error bound). Let t > 0, Aq,. .. ,Ajv S 
and uq G C". Let co{t),... ,Ck-i{t) be the result of Algorithmic let Ek^p{f) 
he defined by ([33]). Then, the total error in the coefficients || vec(£'fc(t))|| sat¬ 
isfies 

II vec(.EA.^p(t))|| = ||exp(tLfc)uo -4(^)1! < 2||no|| 

where 

N N 

(3.6) a = J2\\M\ o-nd fi = p{Ao)+ J2\\M\^ 
e=o e=i 

and p{B) denotes the logarithmic norm defined in (11.51) 

Proof. The result follows directly from [ 8 l Thm. 2.1], and Lemma [ 8 | and 
Corollary llfll in Appendix A. □ 

The coefficient error bound in Lemma U] implies the following bound on 
the error errK,fc,p, via the relation 










Theorem 5 (Krylov error bound). Let errx.fc.p defined in (|3.2p corre¬ 
sponding to applyingp steps of Algorithm\J\to Aq,. .. ,A]\f € anduo G C”. 

Then, 


(3.8) 


errK,fc,p(t, e) < 2 



(^^)Peima^{l,/3} 

-^- hoW 


where a and fi are given in ra . 

Proof. By (13.7p and the Cauchy-Schwarz inequality we have that 


k-l 


errx 


,k,p{t,e) = \\^e\cfit) - Ci{t))\\ 


i=o 



k-l 


\\ce{t) - Ci{t)\[^ 


£=0 


Yec{Ek,p{t))\\. 


The claim follows now from Lemma HI □ 

3.2. A bound on the truncation error. The previous subsection gives 
us an estimate for the error in the coefficient vectors ci{t). To characterize 
the total error of our approach, we now analyze the second term in the error 
splitting (I3T]), i.e., the remainder 


k-l 

errT,fc(t,e) := \\u{t,£) - ^e^C£(t)||. 

i=o 

Lemma M of Appendix gives a bound for the norms of cfit) and can be used 
to derive the following theorem which bounds errT,fc(t, e). 

Theorem 6 (Remainder bound). Let co,ci,. .. be the coefficients of the e- 
expansion (ILI of u{t,e) when N >1. Then, the error errT,fc(t, e) is bounded 
as 

(3.9) errT,fc(t,e) < Ci(t,e) ^ 

1=0 ILAf-l ^ ^r 


where 


(3.10) 


Ci{t,£) 

C2{t,e) 


^|sign(|e|-l) ^t{pi{AQ)+eNa)+C2{t,e)-l 

el'^ eNta. 
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Proof. From Lemma [13] it follows that 


errT,fc(t,e) 




e=k 


l=k 


£=k 


(eNta)^ N~\ 


where C'i(t,e) = QiM^o)+eNa) i||itQ||. 

Setting k = k — N\_^\ and using the bound = (c^)^ < c^^sn{c-i) 

for c > 0, we get 


^ , |£ (eNta)^ nI ^ ^ (|g|'^ejVfa)l'iv1 


— U|sign(|£|-1) 


Af-1 




E E 

j=o e=[±\+j 


(|g|'^ eNtaY 


Using the inequality ED Lemma 4.2] 




£=k 


~kr 


for X > 0, 


the claim follows. □ 

We also give a bound for the special case = 1 since it is in this case 
considerably lower than the one given in Theorem [6] 

Theorem 7 (Remainder bound = 1). Let A^ = 1. Then the remainder 
errT,fc is bounded as 


(3.12) errT,fc(L£)< 


jdM(^o)+|£|||^i||)(||£| ||ij4i ID* 

kl 


Ikoll- 


Proof. From Lemma [ID and (IB.3h we see that cg{t) consists now of one 
integral term which can be bounded by Lemma [12] giving 

Therefore 


en.T,»(M) < E 1^1' IkHDII 


£=k 


£=k 


and the claim follows from the inequality (|3.11l) . □ 
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4. An a posteriori error estimate for the Krylov approximation. 

Although the previous section provides a proof of convergence the final bound 
is not very useful to estimate the error. We therefore also propose the following 
a posteriori error estimates, which appear to work well in the simulations in 
Section [5l 

Let Qpex.p{Hp)ei be the approximation of e^b, ||6|| = 1, by p steps of the 
Arnoldi method. Then, due to the fact that our algorithm is equivalent to the 
standard Arnoldi method, the following expansion holds |21) 

OO 

(4.1) e^6 - Qpexp(iLp)ei = hp+i^p'^ep(pi{Hp)ei A^~^qp+i, 

e=i 


where (pi{z) = (j+e)\ Qp+i (P + 1)®^ basis vector given by the 

Arnoldi iteration. 

We estimate the error of the Arnoldi approximation of exp(tLfc)rIo, i.e., 
the approximation of the vector vec{Ek^p{t)), by using the norm of the first 
two terms in (ilj. This gives us the estimate 

(4 2) « hp+l^p{ep^pl{tHp)elqp+l + ep^p 2 {tHp)el {tLk)qp+i)\\uo\\ 

:= erffc,p(t). 

Then, for the Krylov error errK,fc,p(t,e) in the total error (13.ip . we obtain an 
estimate erfK,fc,p(L s) directly using (|3.7D : 

errK,fc,p(L e) = \\Ek,pit)[l, 

(4.3) = II (4 ® [1, e,..., vec(£;fc,p(t))|| 

Pi 11(4 <8) [l,e,.. errfc_p(t)|| =: erfK,fc,p(t, e)) 

where k = 1 + {N —l)p. Notice that the scalars epipi{tHp)ei and eJip 2 {tHp)ei 
in (14.2D can be obtained with a small extra cost using the fact that [H Thm. 2.1] 


[Ip O] exp 


Hp ei ei1 \ 
0 0 0 
0 1 oj / 


[exp(i4) ipi{Hp)ei + ip 2 {Hp)ei ipi{Hp)ei] . 


Since the a priori bound given by Theorem [6] is rather pessimistic in practice, 
in numerical experiments we only use the Krylov error estimate (14.3D as a total 
error estimate when N > 2. For A = 1, we use also the truncation bound 
given in Theorem [71 i.e., the total estimate is then 


(4.4) 


\u{t,e) -Uk,p{t,e)\\ < errT,A:(t,e) + errK,fc,p(t,e) 

gt4(Ao)+NI|Ai||) (1^1 \\tAi\\)P 


errK,fc,p(L e) + 


P 
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5. Numerical examples. The behavior of the algorithm is now illus¬ 
trated for two test problems: one stemming from spatial discretization of an 
advection-diffusion equation and the other one appearing in the literature m 
corresponding to the discretization of a damped wave equation. 


5.1. Scaling of L^- It turns out that the performance of the algorithm 
can be improved by performing a transformation which scales the coefficient 
matrices. This scaling can be carried out as follows. Let • • • j € 

^nxn corresponding block-Toeplitz matrix of the form ( 12 . 211 . 

Let 7 > 0 and define := diag(l, 7 ,..., 7 ”^“^) 0 /„. Then it clearly holds 

c(t) = exp(tLm) uo = Smexp(tS;;^LmEm) uq 

= Emexp(tLm}uo, 

where is the matrix ( 12 . 2 p corresponding to ^ 0 ) 7 ”^^!) • • • 

Thus, we see that using this scaling strategy corresponds to the changes 

(5.1) e —>■ 76 and Ai — 

when performing the Arnoldi approximation of the product exp(tLm) uq. This 
is also evident from the original ODE (|l.ip . 

The performance of the algorithm appears to improve when we scale the 
norms of coefficients Ai, 1 < < N, such that they are of the order 1 or less. 

To balance the norms, we used the heuristic choice 

(5.2) 7 = max 

This was found to work well in all of our numerical experiments, giving both 
good convergence and a posteriori error estimates. 

We note that scaling has also been exploited for polynomial eigenvalue 
problems, e.g., in [7]. Our scaling (15.2p can be interpreted as a slight variation 
of the scaling proposed in [H Thm. 6.1]. Another related scaling, one for the 
matrix exponential of an augmented matrix, can be found in jT] p. 492]. 


5.2. Advection-diffusion operator. Consider the 1-d advection-diffusion 
equation 


d d 

(5-3) —y{t,x) = a-^y{t,x) + £—y{t,x), 


y(0,x) = yQ{x) 


with Dirichlet boundary conditions on the interval [0,1] and 2 /o(x) = 16 ((1 — 
x)x)‘^. The spatial discretization using central hnite differences gives the or¬ 
dinary differential equation u' = (Aq -|- £Ai)u, u ( 0 ) = uq & M"", where the 
matrices Aq and Ai are of the form 


Ao = 


a 

(Ax)2 


■-2 1 



-1 


1 -2 1 

1 

1 

-1 


1 -2 1 

’ 2Ax 


1 

-1 

1 -2_ 



1 
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where Ax = (n + 1)“^ and uq is the discretization of yo{x). We set n = 200 
and a = 3 • 10“^, and approximate at t = 0.5. Then, ||f^o|| ~ 95. We compute 
the approximations Uk,p{t,s) for e = 10“^, 1.5 • 10“^, and 3 • 10“^. Then, 
respectively, ||te^i|| ~ 0.4, 6.0 and 12.0. Figure [5T] shows the 2-norm errors of 
these approximations and the corresponding a posteriori error estimates using 
(14.4p . We observe superlinear convergence for the error and the estimate. 



Figure 5.1. 2-norm errors of approximations Uk,p{t,e) and the estimate (14.411 when e 
has the values £i = 1 • 10~^, £2 = 1-5 • 10~^ and £3 = 3- 10“^. 


To illustrate the generality of our approach we now consider the case N = 2, 
namely a modification of (|5.3p 

d d‘^ d 

(5.4) —y{t,x) = a-^y{t,x) + £—y{t,x) + eHy{t,l-x), y(0, x) = yo(x); 

the extra term can be interpreted as a non-localized feedback. We set the 
parameter a = 3 • 10““^ and 6 = 2- 10^. The spatial discretization with hnite 
differences gives the ODE u' = (Aq + + £^^ 2 ) u-, 11 ( 0 ) = uo £ where 

Uq and the matrices Aq and Ai are as above, and 


^2 = 6 - 


1 


1 


We compute the approximations Uk,p{t,s) for e = 10“^, 1.5 ■ 10“^ and 3 • 10“^, 
for which respectively, ||te^^ 2 || ~ 5.0 ■ 10“^, 0.11, and 0.45. Figure [5A] shows 
the 2-norm errors of these approximations and the corresponding a posteriori 
error estimates using (14.4p . 
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In Figure [513] we illustrate the dependence of the convergence on the scaling. 
Clearly the choice (15.2p results in the fastest convergence for this example. 
Simulations with 7 larger than what is suggested by (|5.2H did not result in 
substantial improvement of the convergence. 



Figure 5.2. 2-norm errors of approximations Uk,p{t,s) for the equation (15.41) and the 
error estimate gsi) when e has the values ei = 1 • 10 ^, £2 ~ 1-5 • 10 ^ and £3 = 3 • 10 ^. 


5.3. Wave equation. Consider next the damped wave equation inside 
the 3D unit box given in m Section 5.2]. The governing 2?i-dimensional 
hrst-order differential equation is given by 


d 

u{t) 

0 

I 

u{t) 


’«(o)' 


Uo 

dt 


' 1 

1 

-M-1(7(7)_ 

u'{t) 

1 

u'{0) 


u()_ 


where (7(71,72) = 7i(7i + 72(72- The model is obtained by finite differences 
with 15 discretization points in each dimension, i.e., n = 15^. The matrix 
K denotes the discretized Laplacian, (7(71,72) the damping matrix stemming 
from Robin boundary conditions, and M the mass matrix. We carry out 
numerical experiments for parameter values 71 = 0,1, 2 and 72 = 0,1, 2. 

We reformulate (15.5D in the form (II.ip by setting 


Aq = 


0 

-M-^K 


I 


7li = 


0 

-M-1(72 


Then, the variable e in (|l.ll) corresponds to 72. This means that by running 
the algorithm for a fixed value of 71, we may efficiently obtain solutions for 
different values of t and 72. 
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Figure 5.3. 2-norm errors of approximations Uk,p{t,e), when e = 1.5 • 10 ^ using 
different scalings dSIJ. The last option corresponds to the scaling ISM- 


Figures 15.41 show the contour plots of the numerical solutions of (|5.5p at 
t = 9 on the plane {{x, y, z) G [0, 1]^ : z = 0.5} for different values of ( 71 , 72 )- 
Note that for a fixed value of 71 , only one run of the algorithm is required to 
compute the solution for many different 72 . 

In Figure 1^31 we illustrate the relative 2-norm errors of the approximations, 
when 7 i = 2 and 72 =0,1 and 2. Then, ||t^o|| ~ 108, and, respectively, 
||t 72 ^i|| ~ 0,9.6 and 12.9. We again observe superlinear convergence, and, 
moreover, the a posteriori error estimate is very accurate for this example. 

6. Conclusions and outlook. The focus of this paper is an algorithm 
for parameterized linear ODEs, which is shown to have superlinear convergence 
in theory and perform convincingly in several examples. The behavior is con¬ 
sistent with what is expected from an Arnoldi method. Due to the equivalence 
with the Arnoldi method, the algorithm may suffer from the typical disadvan¬ 
tages of the Arnoldi method, for instance, the fact that the computation time 
per iteration increases with the iteration number. The standard approach to 
resolve this issue is by using restarting, which we leave for future work. 
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Appendix A. Technical lemmas for the norm and the field of 
values of Lm- 

We now provide bounds on the norm and the field of values of L^, which 
are needed in Section 13.11 The derivation is done with properties of field of 
values. Recall that the field of values of a matrix A G is defined as 

F{A) = {x*Ax : X G C"', ||x|| = 1}. 

The bounds for the norm and the field of values of An follow from the 
block structure of L^- 
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Lemma 8. Let N >t) and be given by (|2.2D . Then, 

N 

\\Lm\\<Y.\\M\ 

e=o 

Proof. Let x = € C"™' such that Xi G C” for all 1 < i < m 

and ||x|| = 1. From the block Toeplitz structure of L^n we see that 


N 


|^mX|| ^ ^ ^ , 


t=0 


n—i 


N 


\ k=l 


e=o 


\ k=l 


N 


2 — 


Ell All- 


£=0 


□ 

Next, we give a bound for the field of values of the matrix L^- Let d{S, z) 
denote the distance between a closed set S and a point z. 

Lemma 9. Let N > 0 and Lm be given by (12.2p . Then, 

N 

P{Lm) C{zec-. d{F{Ao),z) < Y, Pdl}- 

l=l 


Proof. Let x = [xj' ... where Xj G C” for all 1 < i < m and 

||x|| = 1. Then 

m 

(A.l) X*LmX = X^AqX^ + X*LmX, 
e=i 

where L^ equals Lm on the subdiagonal blocks and is otherwise zero. By the 
convexity of the field of values m Property 1.2.2], we know that YlT=i G 

J^(Ao). The second term in ()A.1I1 can be bounded as in proof of Lemma [8l 
giving 

N 

|x*LmX| < ||Lm|| < Yj 

£=1 


□ 

As a corollary of Lemma [9l we have the following bound, which follows 
directly from the fact that the logarithmic norm of a matrix in 2-norm equals 
the real part of the rightmost point in its field of values. 

Corollary 10. Let Lm be given by (12.2p . Then, 

N 

(A.2) /i(Lm) < ir(Ao) + ||AfII, 

£=1 
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where ^-{A) denotes the logarithmic norm of A defined in (11.51) . 

Appendix B. Technical lemmas for coefficient bounds. We now 

derive bounds needed for the a priori analysis of the truncation error in Sec¬ 
tion 13.21 The following result provides an explicit characterization of the ex¬ 
pansion coefficients. The proof technique is based on the same type of reason¬ 
ing as what is commonly used in the analysis of Magnus series expansions for 
time-dependent ODEs; see, e.g., j^. 

Lemma 11 (Explicit integral form). Let i and N be positive integers such 
that N < i. Denote by Ci the set of compositions of I, i.e., 

(B.l) Cl = {(*1 ,..., ir) G : ii + ir = i}, 

and further denote 

(B.2) Ce^N ■= {(*1 ,... fir) & Cl : is < N for all 1 < s < r}. 

Then, 

co{t) = 

t 

ci{t) = Y. [ f 

(B.3) (b,.'-,*r)eCf,jv 0 0 

—1 

■ ■■ J dtq ... for i > 0. 

0 

Proof From the ODE (12.4p and the variation-of-constants formula it follows 
that 

(B.4a) co{t) = e*^°no, 

min{Ar,£} * 

(B.4b) Ci{t) = ^ ds for £ > 0. 

k=i { 

Using (IB.4p we now prove (IB.3P by induction. For f = 1, we have Ci = {(1)} 
and Ciq = {(!)}■ From (IB.4aP and (IB.4bp we directly conclude that 

ci(t) = [ dti. 

Jo 

Suppose (|B.3P holds for i = 1,... ,p — 1 for some value p > 1. From Defini¬ 
tion we know that the row sum of any element of Cp is p, and the row 

sum of any element in Cp-k is p — k. Therefore, Cp satisfies the recurrence 
relation 

p 

(B.5) Cp=[J IJ {kfii,...fir) 

k 1 (21 , . . 
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and 


min(A^,p) 

(B.6) Cp^N = U U 

k=l {ii,---,ir)&Cp-k,N 

By using (IB.4bp with i = p and the fact that ()B.3h is assumed to be satisfied 
for i = 1,... ,p — 1 we have 

min{Af,p} ^ 

Cpit) = '^2 I ds 

k=i { 

min{Ar,p} * « 

= ... 

k=l Q {ii,---,ir)^Cp-k,N 0 

^ir — l 

j ... dtj^ds. 

0 

By rearranging the terms as a double sum and using (IB.61) . we have 


Cp{t) 


min{7V,p} * ^ 

Y. '22 ... 

k=l {ii,..;ir)&Cp-k,N 0 0 

*V-1 

J ... dti^ds 

0 

t ^V-l 

'22 J ... J 

(n,...,ir)eCp.jv 0 0 


which shows that (IB.311 holds for i = p and completes the proof. 

□ 

Lemma 12. Let m,N be positive integers, N < m, and let Cm,N defined 
as in (|B.2|) . Let {ii,i 2 , ... fir) & Cm, AT, a = maxj=i^,,,_jv ll^jll; assume that 
for all 1 < j < N. Then, one eorresponding term in (lEl is bounded as 


(B.7) 


E ‘'T — L 


0 0 


< e‘ 


r\ 
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Proof. By using the Dahlquist bound [221 p. 138] for the matrix exponential, 
the rightmost integral on the left-hand side of ()B.7I1 can be bounded as 

*v-i 

0 

*V-1 

0 

The claim (|B.7p follows by applying the same bounding technique r — 1 times 
for the remaining integrals, and using that ti <t for any i. □ 

Lemma 13 (Coefficient bound). Let co,ci,. .. be the e-expansion ofu{t,e) 
in (|1.2I) for > 1 m (11.21) . and let a = maxj=i^...^7v ll^jll- Then, for any i >0 
such that k := \^~\ >2, 


Proof We first note that the maximum length of any vector in Cm,N is m, 
and the vector with the shortest length has length at least k = Hence, 

Lemma m can be rephrased as 


m 

Cmit) = 

r=k 


t 


(B.8) 

Lr-1 

0 

Since, Cm,N C Cm we can bound the number of elements in Cm,N 

( TTl — 1 \ 
i 1 ) I 

and Lemma [12] and ()B.8I) imply that 

£=k ^ ^ 

Moreover, 


m — 1 
r — 1 


(m — l)(m — 2 ) • • • {m — r -|- 1) 
(r-1)! 

21 


< 


m 


(r-1)! 










and therefore 


(B.9) 


< e 


tuj{Ao) 


E 

r=k 


rrCitaY 
(r — 1)! r! 


< e' 




E 

r=k 


(r — 1)! r! 


In the second inequality in (IB.9P we use m = = Nk < Nr. 

Using the inequality e < nl, n > 1, we see that for k >2 


r^Nta) 


<e*MAo)y- 


— N^^{Ao)-l 


= e 


(eNtaf Y 


ll-oll < ^ 

r=k 

{eNtaY 




^ (r + A; — 1)! r! 


{eNtaY 

e(r — 1)! 

Ikoll 


Ikoll 


r=0 


t{ii{Ao)+eNa)-lf \k 

< ----^Ikoll, 


(fc-1)! 

where in the last inequality we use < ■^fc:riy!- D 
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